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Abstract 

Reaction networks in the bulk and on surfaces are widespread in physical, chemical and biological 
systems. In macroscopic systems, which include large populations of reactive species, stochastic 
fluctuations are negligible and the reaction rates can be evaluated using rate equations. However, 
many physical systems are partitioned into microscopic domains, where the number of molecules in 
each domain is small and fluctuations are strong. Under these conditions, the simulation of reaction 
networks requires stochastic methods such as direct integration of the master equation. However, 
direct integration of the master equation is infeasible for complex networks, because the number of 
equations proliferates as the number of reactive species increases. Recently, the multiplane method, 
which provides a dramatic reduction in the number of equations, was introduced [A. Lipshtat and 
O. Biham, Phys. Rev. Lett. 93, 170601 (2004)]. The reduction is achieved by breaking the network 
into a set of maximal fully connected sub- networks (maximal cliques). Lower-dimensional master 
equations are constructed for the marginal probability distributions associated with the cliques, 
with suitable couplings between them. In this paper we test the multiplane method and examine its 
applicability. We show that the method is accurate in the limit of small domains, where fluctuations 
are strong. It thus provides an efficient framework for the stochastic simulation of complex reaction 
networks with strong fluctuations, for which rate equations fail and direct integration of the master 
equation is infeasible. The method also applies in the case of large domains, where it converges to 
the rate equation results. 

PACS numbers: 05.10.-a,82.65.+r 



I. INTRODUCTION 



Reaction networks commonly appear in physical, chemical and biological systems, where 
reactions may take place in the bulk or on a surface. When the surface or bulk systems 
are macroscopic, the populations of reactive species are typically large and the law of large 
numbers applies. Thus, fluctuations in the concentrations of the reactants and in the reaction 
rates become negligible. As a result, these reaction networks can be analyzed using rate 
equation models, which account for the average concentrations and ignore fluctuations. 

In some cases, the system is partitioned into small domains, such that the reactants 
cannot diffuse between them. The populations of reactive species in each domain become 
small and their fluctuations cannot be ignored. As a consequence, rate equations fail and 
the simulation of these reactions requires stochastic methods such as direct integration of 
the master equation The master equation takes into account the discrete nature of the 
reactants as well as the fluctuations. It is expressed in terms of the probabilities of having a 
given set of population sizes of the reactive species in a given domain. In certain cases, such 
as radioactive decay, an analytical solution based on generating functions is available 2|. In 
other cases, numerical methods are required. For simple reaction networks that involve few 
reactive species, numerical integration of the master equation is useful and efficient 
However, as the number of reactive species increases, the number of variables in the master 
equation quickly proliferates 6|, |7|, making the direct integraion infeasible. 

Here we focus on networks in which reactions take place between pairs of species, and 
the reaction products may be reactive or non-reactive. Such networks may be described 
by graphs: each reactive species is represented by a node; the reaction between a pair of 
species is represented by an edge that connects the corresponding nodes. Typically, these 
networks are sparse, namely most pairs of species do not react with each other. For such 
sparse networks, the recently introduced multiplane method provides a dramatic reduction 
in the number of equations 8|]. The method is based on breaking the network into a set of 
maximal fully connected subnetworks (maximal cliques). It involves an approximation, in 
which the correlations between pairs of species that react with each other are maintained, 
while the correlations between non-reacting pairs are neglected. The result is a set of lower 
dimensional master equations, one for each clique, with suitable couplings between them. 
For sparse networks, the cliques are typically small and mostly consist of two or three nodes. 
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This method thus enables the simulation of large networks much beyond the point where 
the master equation becomes infeasible. 

The multiplane method has already been used in the simulation of cornplex chemical 



networks on dust grains in interstellar clouds 



where rate equations fail 
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while direct integration of the master equation is impractical [y, |7[ • The multiplane method 
so required for the simulation of genetic networks in cells, where the master equation 
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are not applicable for large networks. 



In this paper we analyze the multiplane method and examine its validity. This is done 
by comparing the results with the complete master equation. The comparison is done both 
for the probability distributions and for the first and second moments, which represent 
the population sizes of reactants and reaction rates, respectively. It is shown that the 
multiplane method provides accurate results for both the population sizes and reaction 
rates. For concreteness, we use below the terminology of surface reactions. In this context, 
the small domains are taken to be dust grains (assumed to be spherical, for simplicity), 
and the reactants are atoms or molecules that enter the system as incoming flux from the 
surrounding gas phase (below we use the words atoms and molecules interchangeably) . The 
reactants and reaction products leave the system by thermal desorption. The reactions that 
take place on a grain are driven by diffusion of reactants on its surface until they encounter 
each other and react. In spite of this specific terminology, the multiplane method can be 
adapted to other contexts, such as reactions in a solution, protein interactions in a living 
cell and birth-death processes in population dynamics. Here we focus on the calculation of 
steady-state solutions and thus do not expand on birth-death processes which may exhibit 
absorbing states. 

The paper is organized as follows. In Sec. [TTl we briefly review the rate equation and 
master equation methods, presenting them for a simple reaction network. In Sec. IIIII we 
describe the multiplane method. In Sec. [IV] we test the performance of the multiplane 
method. An analysis of the method in the limits of small and large grains is presented in 
Sec. |Vl In Sec. |VT] we show how to apply the method to more complex networks. In Sec. 
IVIII we briefly describe its applications in interstellar chemistry and in genetic networks. 
The main findings are summarized and discussed in Sec. IVIIII 
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II. THE RATE EQUATIONS AND THE MASTER EQUATION 



Consider a small grain, exposed to fluxes of three different atomic species, denoted by 
Xi, X2 and X3. The adsorbed atoms on the grain reside in adsorption sites. The number of 
sites, S, is proportional to the surface area of the grain. The incoming flux of Xj, i = 1, 2, 3, 
is given by fi (s~^) atoms per site. Thus, the flux of atoms per grain is Fi = fiS (s~^). 
The adsorbed atoms may desorb due to thermal excitations. The desorption rate of the Xi 
species from the grain is denoted by Wi (s~^). While residing on the grain, the atoms diffuse 
on the surface via hopping between adjacent sites. The hopping rate of X^ atoms is given by 
a, (hops s~^). It is convenient to deflne the sweeping rate Ai = ai/S, which is approximately 
the inverse of the time it takes an adsorbed Xi atom to visit nearly all the adsorption sites 
on the grain sur 



appears in Ref. 
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19|. A more accurate expression for Ai in the case of spherical grains 



201 1 . where it is shown to be reduced by a logarithmic factor. 



The diffusion process enables adsorbed atoms to encounter each other and react. Here 
we consider a simple reaction network that includes the reactions Xi + X2 — > X4 and 
Xi +X3 — > X5, where the X4 and X5 molecules are the reaction products. The graph that 
illustrates this network is shown in Fig. [D^a). 



A. The Rate Equations 

The rate equations that describe the network of Fig. [I](a). take the form 



{A, + A2){N,){N2) - (Ai + As){Ni){Ns) 
{A, + A2){N,){N2) 

{A, + A3){N,){N3), (1) 

where (Ni) is the average population size of X^ atoms on a grain. The flrst terms on the right 
hand side of Eq. ([T]) represent the incoming fluxes of Xi atoms. The second terms represent 
the desorption of Xi atoms, which is proportional to the Xi population on the grain. The 
remaining terms account for the reactions between adsorbed atoms. The production rates 
of X4 and X5 molecules per grain (in units of s~^) are given by R4, = {Ai + A2){Ni){N2) 



d{N,) 

dt 
d{N2) 

dt 
d{Ns) 

dt 



F,-Wi{N,) 

F2-W2{N2) 

F,-Ws{Ns) 
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and i?5 = {Ai + A3){Ni){Ns). For simplicity, we assume that non-reactive product species 
desorb into the gas phase immediately upon formation. 

For large grains, Eqs. ([T]) account correctly for the reaction rates. However, in the limit of 
small grains, some of the average population sizes, (iVj), may become small. In this case the 
discrete nature of the adsorbed atoms and molecules becomes important and the fluctuations 
cannot be ignored. As a result, the reaction rates obtained from the rate equations ([1]) are 
incorrect and stochastic methods are needed. 

We also consider a related network, shown in Fig. [T](b), in which is the product of the 
reaction between Xi and X2 (namely, X3 and X4 are the same species). The rate equations 
that describe this system are the same as in Eq. ([T]) except that in the third equation one 
needs to add the term [Ai + A3){Ni){N3). The production rate of X5 is still given by 
defined above. 



B. The Master Equation 

The dynamical variables of the master equation are the probabilities P(ni, n2, tt-s) of 
having a population of rii atoms of species Xi on the grain. It takes the form 

3 

P(ni,n2,n3) = ^ Fi[P(. . . , - 1, . . .) - P(?t,i, n2, ns)] 

i=l 
3 

i=l 

+ (Ai + A2) [(rii + l)(n2 + l)P(ni + 1, ^2 + 1, rig) - nin2P(ni, n^, n^)] 
+ {Ai + A3)[{ni + l)(n3 + l)P(rai + 1,^2, ris + 1) - nin3P(r2i, ^2, ris)], 

where rii = 0, 1, 2, . . .. The first term in Eq. ([3]) describes the effect of the incoming flux. 
The probability P(. . . , rij, . . .) increases when an Xj atom is adsorbed on a grain that already 
has Hi — 1 adsorbed Xi atoms. This probability decreases when an Xi atom is adsorbed on 
a grain that includes atoms of species Xi. The second term accounts for the desorption 
process. The third and fourth terms describe the reactions that take place on the grain. 

In numerical simulations the master equation is truncated in order to keep the number of 
equations finite. A convenient way to achieve this is to assign upper cutoffs on the population 
sizes of the reactive species, n™'^^, i = 1, . . . , J, where J is the number of reactive species. 
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In the network of Fig. [D^a), the average population size of Xi on a grain is given by the 
first moment 



W= E E E ^^^^''^2,^3). (3) 

711=0 712=0 713=0 

The production rates per grain of and X^ molecules can be obtained from the mixed 
second moments of P(ni, 77-2, na), according to i?4 = {Ai + A2){NiN2) and = {Ai + 
A3){NiN3), where 

max max max 
1 3 3 

{N,N,) = E E E n,n,P{ni,n2,ns). (4) 

711=0 712=0 713=0 

In a network of J reactive species, the number of equations to be solved is 

iVs = n(nr^ + l). (5) 

7=1 

The truncated master equation is valid if the probability to have a population larger than 
the assigned cutoff is negligible. Note that Ne grows exponentially with the number of 
reactive species. This limits the applicability of the master equation to simple networks, 
making it impractical in the case of complex networks which involve many reactive species 

m. 

III. THE MULTIPLANE METHOD 

The recently introduced multiplane method provides a dramatic reduction in the number 
of equations [8|. It thus enables efficient simulations of complex reaction networks. Below 
we describe the method using the network of Fig. [T](a). Note that in this network the 
species Xi participates in both reactions. Since the species X2 and X3 do not react with 
each other, one may assume that for a given population size of Xi, their population sizes 
are almost conditionally independent. Under this assumption, the probability distribution 
of the population sizes can be approximated by [8|] 

P(ni,n2,n3) = P(ni)P(n2|ni)P(n3|ni), (6) 

where P{ni\ni) is the conditional probability that there will be Ui atoms of species Xi given 
that there are ni atoms of species Xi on the grain. 
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In order to derive the multiplane equations, we first insert Eq. into the master Eq. 
([3]), and trace over the population size of X3. Using the fact that Z^ns -^'('^sl'^i) = 1 ^^^^ 
that J2n3 Pi^sl^i) = 0, one obtains 

P(ni,n2) = Fi[P(ni-l,n2)-P(ni,n2)]+F2[P(ni,n2-l)-P(ni,n2)] 
+ Wi[{ni + l)P(ni + 1, 712) - niP{ni, ^2)] 
+ W2[{n2 + l)P(ni, 722 + 1) - niP(ni, ^2)] 

+ {Ai + A2)[{ni + l){n2 + l)P(ni + l,n2 + 1) - nin2P{ni,n2)] 

+ (Ai + A3)[(ni + l)(A^3)ni+iPK + l,n2)-ni(Ar3)„^P(ni,n2)], (7) 

where {N3)n-^ = I]„3 n3P(n3|?7,i). A similar procedure, tracing over the population size of 
X2, leads to the equation 

P(ni,n3) = Pi[P(ni-l,n3)-P(ni,n3)]+P3[P(ni,n3-l)-P(ni,n3)] 

+ Wi[{ni + l)P(ni + 1, ng) - r2iP(r2i, rig)] 

+ Wsliris + l)P(ni, ^3 + 1) - riiP(rii, ris)] 

+ {Ai + A3) [{rii + l)(n3 + l)P(ni + 1, 723 + 1) - riin3P(r2i, ^3)] 

+ (A + A2)[(ni + l)(iV2)„,+iP(ni + 1,773) - 77i(iV2)„,P(77i,773)]. (8) 

These are, in fact, two master equations, one for P{ni, 77,2) and the other for P{ni, ris). These 
two master equations are coupled through the conditional averages {Nj)nj^ where j = 2, 3. 
The conditional average, which is evaluated in each one of these master equations is then 
used, essentially as a rate constant, in the other master equation. The multiplane equations 
are solved by direct numerical integration using standard steppers such as Runge-Kutta. 
At each time step, the probability distributions P{ni,n2) and P(77i,773) are updated. The 
conditional averages {Nj)n-^ are then evaluated and used in the next time step. 

The number of equations is significantly reduced as we replace the three-dimensional set 
of equations for P (771, 712, 773) by two-dimensional sets for P{ni,n2), and P{ni^n^). The 
number of equations in the three-dimensional set is given by Eq. <^ with J = 3. The 
number of equations in each one of the two dimensional sets is (77^"^^ -|- 1)(77™'^^ -|- 1), 7 = 2, 3. 

The multiplane method enables one to calculate the average population sizes {Ni) of all 
the species, as well as the reaction rates, expressed in terms of the second moments (NiNj). 
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Consider, for example, the population size (A^i) of the Xi species. It can be expressed in 
two ways, namely 

(^i)= E EniP(ni,n,) (9) 

ni=0 ni=0 

where i = 2 or 3. In the first case {Ni) is evaluated from P{ni,n2), and in the second case 
it is evaluated from P{ni, ris). A nice property of the multiplane method is that the results 
are identical, as can be seen from Eq. ([6]). The difference is merely in the order in which 
N2 and A^3 are traced out. The multiplane method also provides the reaction rates. For 
example, the production rate of the X4 species [Fig. lH^a)] is given by 

i?4 = {A, + A2)j2 E ^in2P{nu ^2). (10) 

ni=0 712=0 

Note that in the derivation of the multiplane equations, certain dependencies were ne- 
glected. Still, the dependence between all pairs of species that react with each other are 
maintained through the conditional averages, (A^i)ni- These conditional averages are essen- 
tial in order to maintain the desired correlations. If the conditional moments (A'i)^^, in 
the multiplane equations ([7j) and ([H]), are replaced by {Ni) for i = 2,3, these equations are 
reduced, by proper summations, to the rate equations ([1]). In this case all the correlations 
are lost. 



IV. SIMULATIONS AND RESULTS 



To examine the multiplane method we have performed simulations of the reaction net- 
works shown in Fig. [1] The results were compared to those obtained from the complete 
master equation. In Fig. ^a.) we present the average population sizes of the Xi (circles), 
X2 (squares) and X3 (triangles) species on a grain vs. the number of adsorption sites, S, for 
the network of Fig. [I](a), obtained from the multiplane equations under steady state condi- 
tions. In the simulations throughout the paper, we chose to use the parameters Wi = 10~^, 
ai = 10, W2 = 10-3, as = 1, W3 = 10-^ and ag = lO^^ (s^^). This choice reflects the 
mobilities and desorption rates in the network of H, O and OH that appears in interstellar 
grain chemistry 2l|. The production rates of X4 (+) and X5 (x) molecules on a grain, 
vs. S, obtained from the multiplane equations, are shown in Fig. ^h). The results are 
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in excellent agreement with the master equation (solid lines). The rate equations (dashed 
lines) provide accurate results for large grains, but for small grains they show significant 
deviations. We have preformed extensive simulations of this system, using a wide range 
of parameters, and found that the consistency of the multiplane method and the complete 
master equation is always maintained. In the simulations presented above the fluxes were 
Fi = IQ-^S, and F2 = F3 = O.OlFi. 

Note that with the parameters specified above, the incoming flux of Xi atoms is much 
larger than the fluxes of X2 and X^. It is often the case in chemical networks that there 
exists a dominant species, which is more abundant and more reactive than the other species 
(such as hydrogen in interstellar grain chemistry). One could speculate that the dominance 
of Xi is the reason for the remarkable agreement between the multiplane results and the 
master equation results. In order to show that this is not the case, and that the multiplane 
equations are generically applicable, we examine some other parameters. In particular, we 
consider the case in which the flux of Xi is much lower than the fluxes of X2 and X3, namely 
-^2 = -^3 = 10^^ and Fi = O.OI-F2. The population sizes and reaction rates obtained for this 
choice of fluxes are shown in Fig. [3](a) and (b), respectively. Clearly, the excellent agreement 
between the multiplane method and the master equation is maintained in this case as well 
as in all other sets of parameters that we have examined. 

It turns out that the multiplane method applies even when one of the species in one 
clique is a product of a reaction that is included in another clique. To demonstrate this fact 
we consider the network of Fig. [11(b) in which X3 is the product of the reaction between Xi 
and X2. This feature may give rise to some sort of correlation between the population sizes 
of X2 and X3. The question is whether such correlations may reduce the applicability of the 
multiplane method. 

The multiplane equations describing the network of Fig. [11(b) are the same as Eqs. ([7]) 
and (IH]), except that in the last term of the second equation, P(ni + 1,^3) is replaced by 
P{ni + 1, n3 — 1). In FigHl^a) we present the population sizes of the Xi (circles), X2 (triangles) 
and X3 (squares) species on a grain vs. S under steady state conditions, obtained from the 
multiplane method for the network of Fig. [11(b). The production rates of the X3 (+) and X^ 
(x) species are shown in Fig. ID^b) . Even in this case, the multiplane results are in perfect 
agreement with the master equation (solid lines). The rate equations (dashed lines) are 
accurate for large grains but deviate for small grains. Here we chose Fi = 10~^, F2 = O.OlFi 
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and F3 = 0, namely X3 molecules are not accreted from the gas phase, and are produced 
only on the grain. 

Figs. [2] - m demonstrate the usefulness of the multiplane method for the simulation of 
reaction networks on small grains. In particular, it is shown that the multiplane equations 
provide accurate results for the population sizes of reactants, given by the first moments 
(Ni), i = 1,2,3, and for the reaction rates, expressed in terms of the second moments, 
{N1N2) and {N1N3). 

The multiplane method only includes the marginal probability distributions, P{ni,n2) 
and P{ni,n^). However, using Eq. (I6j) one can construct an approximation of the complete 
probability distribution, P(?t,i, 77,2, ns). This approximation takes the form 

PMp(rii,n2,n3) = ^ , (11) 

where the marginal probability distributions P{ni,n2), P{ni,n^) and P{ni) on the right 
hand side are those obtained from the multiplane method. In order to examine the accuracy 
of this approximation, we introduce the deviation distance, 

A= E \Pini,n2,n3) - PMpini,n2,n3)\, (12) 

ni,n2,n3 

which is evaluated under steady-state conditions of the master equation and the multiplane 
equations, where P(ni, n2, tt-s) is the distribution obtained from the master equation. We 
have evaluated A for a range of grain sizes between S = 10^ and S = 10^. It was found 
that in all cases A ^ 1. More explicitly, it varies between A ^ 10~^ to A ^ 10^^ This 
indicates that the reconstructed probability distribution Pmp('^i, ^2, n^) provides a very good 
approximation of P(ni, 77-2, ^3). 

While the second moments which involve pairs of species in the same clique account for 
their reaction rate, such moments for species from different cliques have no direct physical in- 
terpretation. Still, they can be used as an additional test for the accuracy of Pmp('^i, ^2, n^). 
Clearly, one may not expect the multiplane method to provide accurate results for such mo- 
ments because the corresponding correlations are neglected. In Fig. [5](a) we show the 
moment {N2N3) vs. grain size as obtained from the multiplane equations for the reaction 
network of Fig. [T](a) (+). Surprisingly, the results are in agreement with those of the master 
equation (solid line). The corresponding rate equation results for {N2){N3), are also shown 
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(dashed line). In Fig[5](b) we show the moment (A^2^3) vs. grain size, as obtained from the 
multiplane equations (+) for the network of in Fig. [T](b). In this network, the species X3 
is the result of the reaction between Xi and X2, enhancing the correlations between them. 
Indeed, the results of the multiplane method deviate from the master equation results (solid 
line) in the regime of small grains. However, for large grains the results of the multiplane 
method and the master equation coincide and agree with those of the rate equations (dashed 
line). In general, we find that for higher moments of the form (iV"iV2iV|), a,b,c = 1,2,..., 
that involve species from more than one clique, the multiplane method is not reliable in the 
limit of small grains. For large grains the multiplane and master equation results coincide. 

V. ANALYSIS OF THE METHOD 
A. The Limit of Small Grains 

Consider the probability distribution P(ni, n2, ^3) in the limit of small grains, where the 
average population sizes satisfy (Ni) ^ 1 for z = 1, 2, 3. In this limit, (A^^) can be expressed 
by (Ni) ^ Pie, where pj < 1 is a constant that depends on the parameters and e -C 1 is 
proportional to the grain size, S. In this case, P(0, 0, 0) is the highest probability while 
P(l, 0, 0), P(0, 1, 0) and P(0, 0, 1) are of order e. The probability P(0, 1, 1) of having a pair 
of X2 and X3 atoms reside simultaneously on the grain is of order e^. The probabilities 
P(l, 1, 0) and P(l, 0, 1), of having pairs of atoms of species that react with each other reside 
simultaneously on the grain, are reduced due to the reactions and go like e^. Under these 
circumstances, the average population sizes satisfy 



{N2) 



P(l,0,0) + O(e2) 
P(0,l,0) + O(e2) 
P(0,0,l) + O(e2), 



(13) 



while the second moments that determine the reaction rates satisfy 



(iViiVs) ^ P(l,l,0) + C»(e^) 
(iViiVs) - P(l,0,l) + O(e^). 



(14) 
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Using these relations, one can show that in the hmit of small grains, to first order in e, the 
population sizes of X2 and X3 are statistically independent, namely 

P(n2,n3) ~ P{n2)P{n^) + ^(e^). (15) 

To show this relation, one needs to examine three states of -^V^s), namely (n2, n^) = (0, 0), 
(0, 1) and (1, 0). In all other cases, P{n2, n^) goes like the quadratic or a higher degree of e. 
As an example, we show that P{N2 = 0, A^'s = 0) ~ -P(A'2 = 0)P{N^ = 0), to leading order 
in e. To this end, we evaluate the left hand side 

P{N2 = 0, iVa = 0) = P(0, 0, 0) + P(l, 0, 0) + 0{e^), (16) 
and the right hand side 



P(iV2 = 0)P(A^3 = 0) = [1 - P(iV2 = 1) + C(e')] [1 - P(A^3 = 1) + C(e^ 

= P(0, 0, 0) + P(l, 0, 0) + C(e2). (17) 



Clearly, the relation (1151) is satisfied. This result justifies the applicability of Eq. ([6]) in the 
limit of small grains. 

The calculation of mixed second moments for pairs of species that belong to different 
cliques, such as {N2N3) involves probabilities such as P(A^2 = 1, A^s = 1) for states in which 
species that do not react with each other reside simultaneously on the grain. It can be shown 
that these probabilities do not satisfy the relation of Eq. (fTSl) . This result is consistent with 
the fact that the multiplane method does not provide accurate results for this moment, as 
already observed in Fig. [5](b). 

These results further support the conclusion that the multiplane method is suitable for 
the calculation of moments confined to a single clique and is unsuitable for moments that 
combine different cliques. To test this conclusion in detail we define the ratio 

f ^ ^MP('^1,'^2,'^3) /.qX 

r]{ni,n2,n3) = — — , (18) 

P(rai, 722,^3) 

which is equal to 1 where the multiplane method is accurate and deviates from 1 elsewhere. 
In Fig. [6](a) we display the forty highest probabilities, P(ni, 71,2,^3), obtained from the 
master equation for the network of Fig. [I](b) in descending order. The results are for a 
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small grain of 5* = 500 adsorption sites, for which the population sizes of adsorbed species 
are exceedingly small. The probabilities drop off very rapidly, implying that the first and 
second moments are indeed dominated by the few highest probabilities. In Fig. E](b) we 
show the parameter 77, for the same set of probabilities. It is confirmed that the multiplane 
method is valid only for the largest probabilities. Beyond the first few entries, 77 begins to 
fluctuate. In Fig. M^c) we show an enlarged plot of 77, including the first 17 probabilities. In 
this graph the probabilities are labeled. It is found that for those probabilities associated 
with states in which only species from a single clique reside simultaneously on a grain, the 
multiplane method is in excellent agreement with the master equation. For states in which 
species from different cliques reside simultaneously on the grain, significant deviations are 
obtained. The first significant deviation between the multiplane method and the master 
equation is found for the probability P(0, 1, 1), in agreement with the previous analysis. 

The analysis above shows that the multiplane method is valid in the limit of small grains. 
In this limit, the probability distribution is dominated be a few high probabilities associated 
with small population sizes of the reactive species. These dominant probabilities satisfy the 
approximation of Eq. ([6]). Therefore, the population sizes and reaction rates obtained from 
the multiplane method and the master equation are in excellent agreement. 

B. The Limit of Large Grains 

The applicability of the multiplane method in the limit of large grains is not surprising 
because in this limit even the rate equations provide accurate results. As shown above, the 
rate equations can be derrived from the multiplane equations by removing the conditions 
from the conditional averages. The accuracy of the rate equations shows that in the limit 
of large grains the correlations are negligible and the probability distribution P{ni,n2,n^) 
can be factorized into a product of probabilities of single-species. Therefore, the multiplane 
method provides accurate results for any desired moments of the probability distribution. 

VI. THE MULTIPLANE EQUATIONS FOR COMPLEX NETWORKS 

For sparse reaction networks with fluctuations, the multiplane method was found to 
provide a dramatic reduction in the number of equations compared to the master equation. 
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The method provides accurate results for the populations of reactive species and for the 
reaction rates. The method was presented for simple networks which include only three 
species. However, the generalization to more complex networks is straightforward. Consider 
the network shown in Fig. [71 The probability distribution of the population sizes of the 
reactive species in this network is P{ni, . . . ,nj). To derive the multiplane equations one 
first needs to split the network into maximal cliques, or maximal fully-connected subgraphs. 
For the network of Fig. [7] these cliques are: Ci : {Xi,X2}, C2 : {Xi,X3}, C3 : {Xi,X4}, 
C4 : {Xi, X^, Xq} and C5 : {Xi, Xq, Xj}. The next step is to write down the master 
equation for the marginal probability distribution associated with each clique. This can 
be done using either the top-down approach, which is straightforward but tedious, or the 
bottom-up approach. 

In the top-down approach, the master equation for the marginal probability distribution 
associated with a given clique is obtained by tracing over all the species that do not belong 
to this clique. This procedure is repeated for each one of the maximal cliques. 

In the bottom-up approach, the master equation for the internal reactions in each clique 
is constructed first. Then, the coupling terms between cliques, which include the conditional 
averages are added one by one. These terms account for reactions between species, such as 
Xj, which belong to the clique and species, such as Xk, which do not belong to the clique. 
In the master equation, the reaction between Xj and X^. is described by terms of the form 
UjUkPi- . . ,nj,nk, . . .). In the multiplane equation for the given chque, Uk is replaced by 
{Nk)nj and P{. . . ,nj,nk, ■ ■ ■) is replaced by the marginal probability distribution for the 
clique. For example, the resulting equation for the clique Ci is 

P(ni,n2) = Fi[P{ni-l,n2)-P{ni,n2)] 

+ F2[P{ni,n2-l)-P{ni,n2)] 

+ Wi [(ni + l)P(ni + 1, na) - niP{ni, ^2)] 

+ W2[{n2 + l)P(ni,n2 + 1) -n2P{ni,n2)] 

+ {Ai + A2) [{ni + l)(n2 + l)P(ni + 1, ns + 1) - nin2P(ni, ria)] 

+ Ai[{ni + 2){ni + l)P{ni + 2, 712) - ni(ni - l)P(ni, ^2)] (19) 
7 

+ 51(^1 + ^i) [(^1 + 1) {N^)n^+lPinl + 1, 722) - ni{Ni)n,Pini, 7^2)], 

1=3 
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where Eq. (Q is used in order to justify the replacement of (A^i)„i,n2 by {n.i)n^. We find it 
instructive to carry out the procedure for the chque C5 as well. In this clique, the species 
Xi and Xq are both correlated with X5. When tracing over X^ one must maintain both 
correlations, giving rise to the conditional average (A^5)ni,n6- The resulting equation takes 
the form: 



P{ni,nQ,ni) = ^ Fi[P(. . . , - 1, . . .) - P(ni, rie, riy)] 

1=1,6,7 

+ E Wi[{ni + l)P{...,n, + l,...)-n,P{n^,n^,nj)] (20) 

j=l,6,7 

+ {Ai + AQ)[{ni + l)(n6 + l)P{ni + l,n6 + l,n^) - uiUeP (121,12^,717)] 
+ {Ai + A7)[{ni + l)(n7 + l)P(ni + I,n6,n7 + 1) - nir27P(ni, rie, riy)] 
+ (Ae + A7)[{nQ + l)(ra7 + l)P{ni,ne + 1,^7 + 1) - nen7P{ni,ne,n7)] 
+ E Mi^i + 2)(«» + 1)P(- . . , ni + 2, . . .) - ri,(ri, - l)P(ni, ng, nr)] 

i=l,6 
4 

+ E (^1 + + l){N,)n,+iP{ni + 1, ne, ^7) - r2i(Ar,)„,P(ni, rig, nj)] 

1=2 

+ (Ai + A5)[(ni + l)(A^5)„i+i,„gP(ni + I,n6,n7) - ni(iV5)„,,„gP(ni, rig, n7)] 
+ (^5 + AQ)[{ne + l)(iV5).„j,„g+iP(ni,n6 + 1,^7) - n(i{N5)ni^neP{ni,nQ,n'r)]. 

This network has been simulated using both the multiplane method and the complete master 
equation. The results were found to be in excellent agreement Q]. 

VII. APPLICATIONS TO PHYSICAL AND BIOLOGICAL SYSTEMS 

Stochastic simulations are of great importance in a wide range of physical systems. Be- 
low we present two examples of current research areas in which the multiplane method is 
expected to be useful. 



A. Chemical Networks on Interstellar Grains 



The chemistry of interstellar clouds includes gas-phase reactions as well as grain-surface 



reactions 
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23(1 . Due to the microscopic size of the grains, and the low flux, the population 



sizes of reactive species on the grains may be small and exhibit strong fluctuations. Under 



15 



;hese conditions rate equations are not suitable for the simulation of grain-surface chemistry 
9], [lo|, [n, 1^. To account correctly for the reaction rates, stochastic methods such as 
direct integration of the master equation [s], 0, 5| or Monte Carlo simulations o], [lol, 24 1 
are required. The master equation is more suitable for grain chemistry because it consists 
of differential equations, which can be easily coupled to the rate equations of gas phase 
chemistry. Furthermore, the master equation provides the probability distribution from 
which the reaction rates can be evaluated directly, unlike Monte Carlo methods that require 
to accumulate large sets of data and to perform ensemble or temporal averages. For simple 
networks the master equation is efficient and provides accurate results. However, for complex 
networks, the master equation becomes infeasible. In this case, the multiplane method 
provides efficient stochastic simulations. 

Consider the network of Fig. [71 Using the following substitutions Xi H; X2 — 
0H;X3 ^ H3CO;X4 ^ H2CO;X5 ^ HC0;X6 ^ 0;X7 ^ CO, this network coincides 
with the reaction network that leads to methanol production on grains in molecular clouds 



lOj. Current experimental effort is aimed at the evaluation of the relevant rate 



constants for the surface diffusion, reaction and desorption of the species involved in this 



network 
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261]. These experiments include infra-red spectroscopy as well as temperature- 



programmed desorption runs using a mass spectrometer to detect the desorbed molecules. 
The resulting parameters, inserted into the multiplane equations, will enable to evaluate 
the reaction rates in interstellar environments and to compare the results with observations. 
The multiplane method for this network provides a reduction in the number of equations 
from about one million, using the master equation, to about one thousand equations. For 
more complex networks the master equation becomes infeasible while the multiplane method 
remains efficient. 



B. Genetic Networks in Cells 

Another important field in which the multiplane method is expected to be useful is the 
study of genetic regulatory networks in cells. These networks describe the transcription of 
mRNAs from genes and their translation into proteins. The regulation is performed at the 
transcriptional level (by transcription factors that bind to the promoter site of the regulated 
gene), at the post-trancriptional level (e.g., by small non-coding RNAs) and at the post- 
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translational level (e.g., by protein-protein interactions). Analysis of these networks revealed 
modular structure. In particular, modules or motifs which perform specific functions and 



repeatedly appear in different parts of the network were identified 
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28 



291] . Common 



examples of such motifs are the autorepressor 
loop [Mil- Other modules such as the genetic switch 



301 ] and different versions of the feed forward 



321 ] and the mixed-feedback loop 33] 



also appear, but are not as common. 

Genetic networks often exhibit strong fluctuations due to the fact that some of the tran- 
scription factors appear in low copy numbers. Moreover, the transcriptional regulation 
is typically performed by a small number of transcription factors which bind and unbind 
to the promoter site at a fast rate. This gives rise to strong fluctuations in the tran- 
scription rate of the regulated gene. Some modules, such as the autorepressor, the ge- 
netic switch and the mixed-feedback loop include positive or negative feedback mechanisms, 
which tend to enhance the role of fluctuations. In particular, the dynamics of the genetic 
switch system was studied extensively using both deterministic and stochastic methods 



34 



35 



37 



38 



39 



40 



41| . It was found that fluctuations play a crucial role in this 



system. While the analysis of small modules such as the genetic switch can be done using 
the master equation, it quickly becomes infeasible when larger networks are considered. The 
implementation of the multiplane method in this context is expected to provide a broader 
perspective on the role of fluctuations in genetic networks. Recently, such fluctuations at the 
level of single cells were measured experimentally using the green fluorescent protein 42]. 
Such measurements are also expected to provide the effective rate constants of the relevant 
processes in the cell. 



VIII. SUMMARY AND DISCUSSION 



We have shown that the multiplane method provides efficient simulations of complex re- 
action networks with fluctuations, for which the rate equations fail and the master equation 
is infeasible. The multiplane equations are obtained by breaking the network into maximal 
cliques and writing down the set of master equations for the marginal probability distribu- 
tions of these cliques, with a suitable coupling between them. For typical sparse networks, 
the method provides a dramatic reduction in the number of equations. We found that the 
multiplane results for the first and second moments, which account for population sizes and 
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reaction rates, respectively, are in excellent agreement with those of the complete master 
equation. It also accounts correctly for higher moments which involve species from the same 
clique. However, the method does not account correctly for second and higher moments 
which include species from different cliques. 

The numerical results are complemented by an asymptotic analysis of the small and large 
grain limits. A more rigorous analysis shows that the multiplane method is asymptotically 



exact in both limits 



43| . It is performed in a more general setting, in which the maximal 



cliques may be broken into smaller cliques. In particular, one may break the entire network 
into cliques of two species each. It is shown that even in this case, in the limits of small and 
large grains, the method still provides exact results for all the first moments and for those 
second moments that involve species in the same clique. 

A related approach, based on moment equations, also provides efficient stochastic sim- 
ulations of reaction networks 4^. In this approach, one constructs differential equations 



for the first and second moments of the probability distribution. The number of equations 
is further reduced to one equation for each reactive species (node) and one equation for 
each reaction (edge). Thus, for typical sparse networks the complexity of the stochastic 
simulation becomes comparable to that of the rate equations. In applications such as inter- 
stellar chemistry, in which the main objective is to calculate the reaction rates, the moment 
equations are advantageous. However, in systems such as genetic networks, in which the 
probability distribution itself is of interest, the multiplane method is required. 
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FIG. 1: Graphic representations of two reaction networks that involve three reactive species. The 
nodes represent reactive species and the edges represent reactions between pairs of species. The 
reaction products are specified near the edges. In these networks there are two cHques: one consists 
of Xi and X2 and the other consists of Xi and X-^. (a) The reaction products, and X^ are 
non-reactive; (b) The product of the reaction between Xi and X2 is the reactive specie X3. 
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FIG. 2: (a) The population sizes of the Xi (circles), X2 (triangles) and X3 (squares) species per 
grain vs. the number of adsorption sites, S, on the grain, obtained from the multiplane equations, 
for the network shown in Fig. [TJa). The results are in perfect agreement with the master equation 
(solid lines) and the rate equations (dashed lines); (b) The production rates of (+) and X^ 
(x) molecules per grain vs. 5, obtained from the multiplane equations. The results are in perfect 
agreement with the master equation (solid lines). For small grains, the rate equation results (dashed 
lines) for the reaction rates exhibit large deviations. Here, Xi is the dominant species. 
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FIG. 3: The population sizes (a) and the production rates (b) per grain vs. S for the network 
of Fig. dja). The multiplane results (symbols) are in perfect agreement with the master equation 
(solid lines). The rate equation results (dashed lines) deviate significantly for small grains. Here 
the species Xi is dominated by X2 and X3. 
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FIG. 4: The population sizes (a) and the production rates (b) per grain vs. S, for the network of 
Fig. in which is the reaction product of Xi and X2. The multiplane results (symbols) are 
in perfect agreement with the master equation (solid lines). The rate equations (dashed lines) are 
also shown. 
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FIG. 5: (a) The second moment {N2N3) of P{ni,n2,ns) vs. S, obtained from the multiplane 
method (+) for the network shown in Fig. [H^a). This moment is not related to any reaction rate, 
thus the multiplane method is not designed to approximate it well. Still, it turns out that the 
results are in good agreement with the master equation (solid line). The rate equation results for 
the corresponding term, {N2){N-i), are also shown (dashed line); (b) The moment, {N2N3), vs. S 
for the network shown in Fig. mb). The multiplane results (+) deviate from those of the master 
equation (solid line) in the limit of small grains. For large grains the multiplane results coincide 
with the master equation and with the corresponding term, {N2){N3), of the rate equations (dashed 
line) . 
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FIG. 6: (a) The probabilities P(ni, 712,^3), obtained from the master equation, arranged in 
descending order from the largest (left) to the smallest (right) for a small grain with S = 500 
adsorption sites. Here we show the forty largest probabilities; (b) The ratio parameter, r], defined in 
Eq. (jlSp . between the probabilities obtained from the multiplane equations and the corresponding 
probabilities obtained from the master equation. The ratio is unity for the first few probabilities 
and then it fluctuates for the the rest; (c) An enlargement of the first seventeen probabilities 
displayed in (b). Standing out are -P(0, 1, 1) and -P(l, 1, 1), for which the multiplane equations and 
the master equation differ. These probabilities have no significant effect on the production rates 
and population sizes, but are expressed when computing other moments. 
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FIG. 7: A graph that represents a complex reaction network which consists of seven reactive 
species. In the multiplane method, this network is broken into five cliques, and a lower dimensional 
master equation is constructed for the marginal probability distribution associated with each clique. 
In interstellar chemistry, this is the network leading to the formation of methanol on dust-grain 
surfaces. 
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